Delocalization induced by nonlinearity in systems with disorder 
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We study numerically the effects of nonlinearity on the Anderson localization in lattices with 
disorder in one and two dimensions. The obtained results show that at moderate strength of 
nonlinearity a spreading over the lattice in time takes place with an algebraic growth of number of 
populated sites An oc f". This spreading continues up to a maximal dimensionless time scale t — 10^ 
reached in the numerical simulations. The numerical values of v are found to be approximately 
0.15 — 0.2 and 0.25 for the dimension d = 1 and 2 respectively being in a satisfactory agreement 
with the theoretical value d/{3d + 2). During the computational times t < 10^ the localization is 
preserved below a certain critical value of nonlinearity. We also discuss the properties of the fidelity 
decay induced by a perturbation of nonlinear field. 

PACS numbers: 05.45.-a, 03.75.Kk, 05.30.Jp, 63.50.-x 
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I. INTRODUCTION 

The phenomenon of the Anderson localization [l| in 
systems with disorder has been extensively studied for 
electron transport and linear waves (see e.g. 0). A re- 
markable experimental progress with the Bose-Einstein 
condensates (BEG) in optical lattices (see e.g. reviews 
S 0] ) stimulated the interest to investigations of the 
effects of nonlinearity on localization. At present the sig- 
natures of localization of BEG in one-dimensional (ID) 
optical disordered lattices have been detected by different 
experimental groups 0, H, M, [13 ■ The effects of non- 
linearity appear also for experiments with BEG in kicked 
optical lattices [111, [H, [3 where the quantum chaos in 
the Ghirikov standard map (kicked rotator) is inves- 
tigated. A similar type of problem also comes out for 
propagation of nonlinear waves in disordered photonic 
lattices which are now actively studied experimentally 
[H, [ll]. In addition to that the problem of lasing in 
random media is also linked to the interplay of lo- 
calization and nonlinearity that makes it related to the 
important field of nonlinear wave propagation in disor- 
dered media [3 ■ In this work we concentrate our studies 
on the time dependent wave packet spreading in presence 
of disorder and nonlinearity leaving aside the problem of 
directed flow and scattering in nonlinear media (see e.g. 
Rcfs. in [l3 and more recent p^). 

The theoretical treatment of the interplay between lo- 
calization and nonlinearity uses numerical simulations 
(see e.g. 0, M M MMMM) and 
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various ana- 



3(1 V However, even 



lytical tools (see e.g. ^[^^mi,^.^ , 

rather powerful analytical tools "Mil, an do not 
allow to obtain the full solution of this rather complex 
problem. The existing rigorous mathematical results 
show that for a sufficiently small nonlinearity there exists 
a Kolmogorov-Arnold-Moser intcgrable localized regime 
for almost all initial conditions [311 but these results are 



" jhttp : //www ■ quantware ■ ups-tlse . f r/dlma | 



applicable only to unrealistically small strength of nonlin- 
earity. Due to that the numerical simulations become es- 
pecially important for investigation of this problem. For 
numerical studies it is especially convenient to use a dis- 
crete lattice that allows to push numerical simulations 
to extremely large times. In addition to that the time 
evolution on a lattice is closely linked to the problem of 
energy propagation in complex molecules, e.g. proteins, 
where nonlinear couplin gs give transitions between lo- 
calized linear modes [2J, [2^, [S^l ■ One of the examples of 
such a nonlinear oscillator chain is the Frenkcl-Kontorova 
model in the pinned phase where the linear sound modes 
are localized in space [33l |. 

Recently, the interplay of the Anderson localization 
and nonlinearity has been investigated by a number of 
mathematical methods where a certain number of in- 
teresting mathematical results has been obtained [13, 
m, [3^. However, these methods still should be devel- 
oped further to understand the asymptotic properties of 
spreading in the lattice at moderate strength of nonlin- 
earity. 

In this paper we further develop the old [13, HH and 
recent studies [2^ of the discrete Anderson nonlinear 
Schrodinger equation (DANSE) and present large scale 
numerical simulations of this model in one and two di- 
mensions d (ID, 2D). In addition we perform numerical 
simulations for the kicked nonlinear rotator model (KNR) 
introduced in [2l| . Our numerical results obtained on di- 
mensionless time scales up to < = 10® show that at mod- 
erate nonlinearity, above a certain threshold, the wave 
packet spreads unlimitedly over the lattice in such a way 
that the squared displacement of the packet on the lat- 
tice grows according to the algebraic law oc with 
the exponent a w 0.3 — 0.4 for d = 1 and a « 0.25 for 
d = 2. This dependence is in a satisfactory agreement 
with the ID estimates [2l| which give a = 2/5 and the 
analytical estimates of this paper which give a = 1/4 for 
2D. We also study the fidelity decay which shows inter- 
esting properties for the nonlinear evolution described by 
our model. 

The paper is composed as follows: in Section II we give 
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the model description and present simple estimates; the 
results for ID and 2D are presented in Sections III and 
IV respectively, the properties of nonlinear fidelity decay 
are discussed in Section IV; the conclusions are given in 
Section V. 



II. MODEL DESCRIPTION AND ANALYTICAL 
ESTIMATES 

Our system is described by the DANSE model: 

'^^-Qf = Eni'n + /?| V'n I V'n, + V {^n+l + 'ipn-l) , (1) 

where P characterizes nonlinearity, y is a hopping matrix 
element on nearby sites, on-site energies are randomly 
and homogeneous distributed in the range —W/2 < En < 
W/2, and the total probability is normalized to unity 
I 1^ = 1- Here, n is the lattice index, in ID it is 
an integer, in 2D it is an integer vector of lattice indexes 
n = {rix, Uy). For /3 = and weak disorder all eigenstates 
are exponentially localized with the localization length 
I « 96(F/VF)2 (ID) at the center of the energy band and 
In^ ^ iy/WY in 2D [13] ■ Hereafter we set for conve- 
nience h = V = 1, thus the energy coincides with the 
frequency. We emphasize here that the DANSE ([T|) ex- 
actly describes recent experiments with one-dimensional 
disordered waveguide lattices (cf. Eq. (1) in [iBl), and 
it also serves as a paradigmatic model for a wide class 
of physical problems where interplay of nonlinearity and 
disorder is important. The DANSE can be considered as 
the Gross-Pitaevskii equation (GPE) Q taken on a dis- 
cretized lattice. In ID this model was studied recently in 

MM- 

To understand the evolution properties of system ([1]) 
it is convenient to expand tpn in the basis of localized 
eigenmodes at = [2l|: A„ = Y^m Qn,mCm where A 
and C are amplitudes in the basis of sites and eigen- 
modes respectively. Due to the localization of linear 
eigenmodes with length / we have for the transformation 
matrix Qnm exp(-|n - m\/l - ixn.m), where x 

are some random phases. From ^ it follows that the 
amplitudes C in the linear cigenbasis are described by 
the equation 

where Cm are the eigenmode energies. The transitions 
between linear eigenmodes appear only due to the non- 
linear /3-term and the transition matrix elements are 

[2l] |. There are about l^"^ random terms in the sum in 
© with V - r^''/^ so that we have idC/dt - PC^. We 
assume that the probability is distributed over An > l"^ 
states of the lattice basis. Then from the normalization 
condition we have Cm l/(An)^''^ and the transition 
rate to new non-populated states in the basis m is F 



/3^[C|^ ^ /3^/(An)^. Due to localization these transitions 
take place on a size I and hence the diffusion rate in the 
distance AR (An)^/'' of d~ dimensional m— space 
is d[ARf/dt - - /32/2/(An)3 ^ pH^/iARf^. At 
large time scales AR ^ R and we obtain 

An^R"^ (^;)2d/(3d+2)^d/(3d+2) 

Thus as in [2lj for d = 1 we have R^ cx t^/^ and for 
d = 2 this gives R^ oc t^/^^ while for large d the scaling 
is independent of d: An cx t^/'^ . 

The relation ([3]) assumes that the dynamics of nonlin- 
ear chain ([T|) is chaotic. On a first glance it seems that 
it cannot be the case since as soon as An grows with 
time the nonlinear frequency shift 5u) ^ /?|'0nP P/ An 
decreases. However, the physical importance relays not 
on the shift value itself but on its ratio to a frequency 
spacing between frequencies of excited modes which is 
Auj ^ 1/An. The latter relation results from the fact 
that all the frequencies are distributed in a finite en- 
ergy (frequency) band and therefore An states excited 
inside such a band have energy and frequency spacing 
Auj ~ 1/An. The dynamics is chaotic if the overlap pa- 
rameter S = 5u} / Alo ^ P > Pc ^ I. It is important to 
stress that S is independent of An. By its nature this 
criterion is somehow different from the usual Chirikov 
resonance-overlap criterion [3] since in our case the un- 
perturbed system is represented by a set of linear oscilla- 
tors while in [3] the oscillators are nonlinear. However, 
the condition Suj > Auj looks rather natural since in the 
opposite limit Suj <C Auj the coupling between modes 
is very weak if the linear frequencies are linearly inde- 
pendent (that should be true in a disordered potential). 
In addition the investigations of three nonlinear oscilla- 
tors with nonlinear couplings performed in [40| indeed 
confirmed the criterion 5* > 1. Therefore from the crite- 
rion S > 1 we obtain that above some critical nonlinear 
coupling p > Pc ^ const the dynamics remains chaotic 
even if probability spreads over larger and larger parts 
of the lattice. This spreading should follow the relation 
(l3|). During this process the local Lyapunov exponent 
X 5uj ^ PI An decreases to zero since the system size 
is unlimited but locally the dynamics is chaotic. 

Another argument in favor of unlimited spreading can 
be obtained on the basis of certain similarities and par- 
allels with the Frenkel-Kontorova chain. In this non- 
linear chain the number of configurations static in time 
{di}jn/dt = in ([1])) grows exponentially with the length 
of the chain while the energy splitting between these con- 
figurations drops exponentially with the chain length (see 
e.g. (ssj). Therefore, due to this energy quasi-degeneracy 
between these static configurations, it is rather natural 
to expect that during the time evolution a spreading over 
all these configurations continues unlimitedly. 

For P > Pc this spreading corresponds to a regime 
of strong chaos with mixing of all modes. The situ- 
ation for P < Pc may have other mechanisms of slow 
chaos with slower spreading and should be analyzed sep- 
arately. For example, the typical spacing in the resonant 
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terms in Eq. ^ is A2 ^ Em + -B^^ - - ^ 
l/P*^ and it is smaller than the coupling matrix element 
f3Vmm,m,m, ^ P/l"""'^ for pi'''^ > 1. Therefore, it is 
possible that for < /? < /3c ~ 1 there may be 

a propagation of two-modes-pairs on a distance much 
larger than / in a certain similarity with a quantum dy- 
namics of the two interactin g p articles (TIP) in a ran- 
dom potential discussed in [33 • Indeed, Eq. ^ can 
be viewed as a mean field approximation for the TIP 
Hamiltonian considered in [38[. In analogy with the 
TIP problem it is possible to expect that the distribu- 
tion of the probability in the basis of linear eigenmodes 
Cm will be characterized by the Breit-Wigner shape: 

Wm,m, = \Cm{t)CmAt]\Lr- ^/[{E - ^rn - CmJ^ + rV4] 

(see e.g. discussion in [33 for TIP). Here, the value of F 
is given by the above estimates. However, the verification 
of this Brcit-Wigner relation requires further numerical 
tests with a projection on the cigcnbasis of linear modes 
that was not done in this work. In this paper we concen- 
trate our studies on the regime /? ~ 1 > /3c- 



III. NUMERICAL RESULTS IN ID 



momentum representation and in the each representation 
the integration is performed exactly up to the computer 
double precision. 
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FIG. 1: (Color online) Time dependence of the averaged sec- 
ond moment (An)^ in ID for the disorder strength W = A aX 
/3 = (black) and /3 = 1 (red/gray). The straight line shows 
the fit for 100 < t < 10** with the slope en = 0.325 ± 0.003. 
Here and below the logarithms are decimal. 



To test the above theoretical predictions we perform 
the numerical simulations of the time evolution given by 
Eq. ([1]). The split operator scheme on a time step At is 
used for integration: 

{t + At) = QVi) it) , (4) 

where Q = exp{—i{En,m + P\ V'n,m \'^)At) is a diagonal 
operator in the lattice space and the application of the 
hopping operator V is done by the fast Fourier transform 
to the conjugated space where V becomes diagonal taking 
in 2D the form V = exp(— 2iAt(cos -I- cos6ny))- For 
the results presented in next Sections we used At = 0.1 
and the averaging was done over Nd = 10 realisations of 
disorder. We checked that a variation of At by a factor 
2-4 does not affect the numerical data for short times 
(e.g. t < 20) and on the large time scales t ~ 10^ the 
statistical behavior of the results remains unchanged. Of 
course, on large times the exact values of tpn are dif- 
ferent for different At due to exponential instability of 
dynamics. But the integration scheme (|4]) is symplec- 
tic and preserves the total probability exactly while the 
total energy is preserved approximately with the accu- 
racy of 1%. Indeed, the final integration step generates 
high frequency ujint = 27r/At « 60 that is significantly 
larger than the energy band width of the linear problem. 
We checked that for ID this integration scheme gives the 
same results as other schemes used in [2^. The lattice 
size in ID was N = 2" site and in 2D iV = 256 x 256. 
The initial state was chosen with all probability on one 
site in the middle of the lattice. We note that for the 
KNR model (O the integration precision is on the level 
of double precision of the computer since the integration 
is done by the fast Fourier transform from coordinate to 



2 p 
1.5 - 



0.5 - 

I ' ' ' ' ' ' ' 

012345678 

log t 



FIG. 2: (Color online) Time dependence of the averaged IPR 
^ for the parameters of Fig. [T]for /? = 1 (red/gray curve) and 
(3 = (black curve). The straight line shows the fit with the 
slope V = 0.125 ±0.001. 

To characterize the properties of time evolution we 
compute one-site probability Wn = IV'tiP, second mo- 
ment of the probability distribution (An)^ in ID and 
(Ai?)^ = (Ahj;)^ + (Ariy)^ in 2D, and the inverse partic- 
ipation ratio (IPR) C = l/Sn''^"'^ which gives an effec- 
tive number of sites populated by the wave packet if all 
Wn^ probabilities arc of the same order. To suppress fluc- 
tuations the quantities (An)^, (AR)^, ^ — l/X^n^"^ 
are averaged over time intervals which are equally spaced 
in logt. In addition the logarithms of these quantities 
are averaged over N^, = 10 disorder realizations. The de- 
pendence on time is fitted by the algebraic dependencies 
(An)2 ~ pi, (Ai?)2 ~ ,c ^ with the exponents 
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The numerical results presented in Figs. ll|2l give values 
of exponents ai = 0.325 ± 0.003 and u = 0.125 ± 0.001. 
The value of ai is in agreement with the data obtained 
in [2^ where it was found that ai = 0.306 ± 0.002 at 
W = 4. Indeed, it should be noted that the standard 
deviation Aai « 0.014 for fluctuations in the exponent 
ai from one disorder realization to another (25j is larger 
than the the formal standard error of the fit of averaged 
data (with Aai p» 0.003). If all states inside the width 
An are populated in a homogeneous way then we should 
have (An)^ ~ and ai = 2i/. The numerical data of 
Figs. I1I2I give this ratio to be aijlv = 1.3 instead of 1. 
This indicates that the probability inside the width An is 
distributed in inhomogcneous way. It is possible that the 
spreading has certain multi-fractal properties that give 
deviations from usual relations between high moments. 
At the same time our data clearly confirm that the IPR 
grows in unlimited way with time (see Fig. [5]). This is 
different from the claim presented in [23|. We attribute 
this difference to the fact that in [23] the data have been 
presented only for one disorder realization and no data 
for fits and their statistical accuracy have been given. At 
the same time the data of [l^] for the second moment 
(An)f are consistent with the results presented here and 
in [111. 

To obtain more results for larger times we also per- 
formed numerical simulations for the KNR model intro- 
duced in [21I . It time evolution is described by the map 
for the wave function: 

V'„(t + 1) = e-*^"'/2-*/5|'/'"l'e-»'=^°'' V„(t) , (5) 

where (n, Q) are the conjugated operators with the com- 
mutation relation [n, G\ = —i and -0 is periodic in B. 
For /3 = this is the model of kicked rotator where all 
quasienergy eigenstates arc ex pon entially localized with 
a localization length I w fc^/2 jl4l. [4l]|. The propagation 
operator is similar to the one of HI) with At = 1, due 
to that it is possible to perform t = 10^ map iterations 
of ^ for the same CPU time as for (jl]). The numerical 
results arc presented in Figs. 3-5. 

These results show that unlimited spreading of proba- 
bility over the sites n takes place at moderate values of 
/3 ~ 1. The probability distribution over n has a plateau 
followed by exponential tails, inside the plateau the prob- 
ability is homogeneously distributed and the width of the 
plateau grows with time (see Fig. [3]). The second mo- 
ment of the distribution and the IPR grow algebraically 
in time with the exponents a\ = 0.387 ± 0.003 and 
V = 0.210 ±0.002 respectively (Figs. IHl). In view of sta- 
tistical fluctuations we consider that these values are in 
a good agreement with the theory estimates (see Eq. ^ 
and [12]). The relation a\ = 2v also works with a rel- 
atively weak deviation from the theory. For the KNR 
model the agreement with the theory is better than for 
the model ([T|). The possible reason is that in the KNR 
all linear eigenmodcs have the same localization length 
I K, fc^/2 while for the DANSE the localization length 
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FIG. 3; (Color online) KNR data ^ for probability distri- 
bution Wn at times t — 10^ (red/gray); l(f (green/gray), 10^ 
(blue black) for P — (top panel with overlapped curves) 
with basis size N = 2^°, for P = 0.03 (middle panel) with the 
same basis size and for 13 = 1 (bottom panel, curves are from 
bottom to top at \n\ x 200 for time from t = 10^ to 10^) at 
iV = 2"; here fc = 3, T = 2. 



depends on the energy value inside the energy band that 
gives stronger statistical fluctuations and require longer 
times for the observation of the asymptotic algebraic 
growth. Also, it is possible that the stronger deviations of 
the exponent values from the theory in ID DANSE model 
are related to the absence of good diffusive approxima- 
tion for the ID Anderson model while for the KNR model 
the diffusive approximation works rather well. 

At small values of /3 = 0.03 the probability distribution 
remains localized during enormously long times t < 10® 
but for larger time t ~ 10^ the distribution grows slightly, 
also ^ and (An)^ are increased by a factor 2 and 3 re- 
spectively (see Figs. 131415^ . It remains unclear if this is 
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FIG. 4: (Color online) The spreading of the second moment 
(An)^ with number of map iterations t for the KNR model 
for the parameters of Fig. O /3 = (black curve), A'' = 2^^°; 
f3 = 0.03 (dashed green/light gray curve), A'^ — 2^'' and f3 = 1 
(red/gray curve), A'^ — 2^^. The straight line shows the fit for 
100 < t < 10^ with the slope ai = 0.387 ± 0.003. 
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FIG. 5: (Color online) Same as in Fig. [5] for the IPR ^ in the 
KNR model. The straight line shows the fit for 100 < f < 10^ 
with the slope u = 0.210 ± 0.002. 

a fluctuation or if there is a very slow (logarithmic ?) 
spreading. In any case the behavior for small nonlinear- 
ity (3 < Pc ^ 0.03 is qualitatively diflFerent compared to 
the case of moderate values oi (3 ^ \ being in a qualita- 
tive agreement with the theoretical expectations. 



IV. NUMERICAL RESULTS IN 2D 

Here we present results for the model ([T]) in 2D. All 
results are averaged over Nd = 10 disorder realisations. 
The time evolution of the probability distribution Wn^.uy 
is shown in Fig. [6] for = 10. At /? = the proba- 
bility is localized while at /3 = 1 it slowly spreads over 
the lattice. The second moment of the space displace- 
ment (Ai?)^ = (Aux)'^ + (Ariy)^ as a function of time 
is shown in Fig. [T] The growth is well described by the 
algebraic dependence (Ai?)^ — Dt"^ . The fit gives the 
values of the exponent 012 = 0.236±0.003 for W = 10 and 
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FIG. 6: (Color online) Probability distribution Wn^,ny inside 
the square 256 x 256 at W = 10 for /3 = (left column) 
and 13 = 1 (right column) and time t = 10^ (bottom pan- 
els), 10® (middle panels); probability is proportional to color 
with maximum at red/gray and zero at blue/black. Top pan- 
els show the decimal logarithm of the integrated probabil- 
ity = E„^W"x,n« for -128 < < 127 at t = 10* 
(red/gray) and 10® (black). 

0.229 ± 0.003 for W = 15. Taking into account that the 
growth of (Ai?)^ is rather slow and that there arc fluctu- 
ations related to disorder averaging the agreement of the 
exponent 02 with the theoretical value a = 1/4 ^ can 
be considered as rather good. In addition the value of 02 
in 2D is decreased compared to the value ai in ID. Their 
ratio a2/ai — 0.233/0.325 = 0.717 is rather close to the 
theoretical value 5/8 given by Eq.(I3]). According to @ 
the ratio D{W = 10)/D{W = 15) = (^(1^ = 10)/1{W 
15))i/2 = {{An{W = 10))l/{An{W = 15))§))i/4 « 1.9 
where (An)g are the values taken at /5 = 0. From the 
data of Fig. [7| at /3 = we have this ratio to be 1.9 
while from data at /3 = 1 we obtain its value as 5 that 
can be considered as satisfactory taking into account all 
fluctuations. 

The time dependence of the IPR ^ is shown in Fig. [S] 
The flt gives the algebraic growth with the exponents 

= 0.282 ± 0.002 for 1^ = 10 and 1/ = 0.247 ± 0.005 for 
= 15 that is in a good agreement with the theoretical 
value 1/4 ([3]) . We note that in 2D the exponents a2 and 

become rather close. This indicates that multi-fractal 
effects become less pronounced in 2D. 

Finally, in Fig. [3] we present the comparison of the 
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FIG. 7: (Color online) Average of the squared spreading 
(AR)'^ as a function of time for two different values of the 
disorder strength W — 10 (top red/gray and black curves) 
and W — 15 (bottom red/gray and black curves) for /3 = 1 
(red/gray curves) and (3 = (black curves). The slopes of the 
straight hue fits for 100 < t < lO" give Q2 = 0.236 ± 0.003 for 
W = 10 and 0.229 ± 0.003 for W = 15. The lattice size is as 
in Fig. H 
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FIG. 8: (Color online) Same as in Fig. [7] but for the IPR C 
The slopes of the straight line fit are v = 0.282 ± 0.002 for 
= 10 (top red/gray and black curves) and f = 0.247±0.005 
for W = 15 (bottom red/gray and black curves); black curves 
are for (3 = 0, red/gray curves are for /3 = 1. 



behaviors of linear system /3 ~ and the one at weak 
nonlinearity /3 = 0.033. These data show that for (3 < (3c 
the behaviors of the two systems are rather similar, for 
times explored in our numerical simulations, that is in 
agreement with the theoretical expectations described in 
Section II. According to our data f3c > 0.033 in 2D. 

To characterize the nonlinear evolution of system ([1]) 
in an additional way we computed another characteris- 
tics which we call nonlinear fidelity defined as f{t) = 
\{'4'n,m{t)\ipn,m}it)\'^ , whcrc ipn,m{t) IS & Small perturba- 
tion of V'n,™ at t = 0. For the linear system with /3 = 
the fidelity f{t) remains constant during time evolution. 
However, for nonlinear dynamics f{t) starts to depend 
on time. Indeed, the perturbation changes the nonlinear 
potential and the system starts to evolve with a slightly 
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FIG. 9: (Color online) Average of the squared spreading 
(A_R)^ as a function of time for for /3 = (black curve) and 
P = 0.033 (red/gray curve) at W = 15. 

different effective Hamiltonian that leads to a decrease 
of fidelity. Such a behavior reminds the fidelity decay 
studied in systems of quantum chaos (see e.g. review 
[4^1). We note that rece ntly the fidelity decay in the 
GPE has been studied in [43], however, there the ampli- 
tude of random potential was considered as a very small 
perturbation while in our case the disordered potential 
is strong and plays a dominant role. Also in [43| the fi- 
delity was considered for perturbation of potential while 
we consider the perturbation of nonlinear field that was 
not addressed in [i^. 
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FIG. 10: (Color online) Dependence of averaged /(f)//max as 
a function of the rescaled time Ft for W = 8, /3 = 1 and 7 
different values of perturbed probability 0.001 < SP < 0.03. 
The average is done over Nd = 12 disorder realizations. The 
straight line shows the dependence /(t)//max = exp{—Tt). 

To study the properties of f{t) we start at i = from 
one lattice state for -0 while the state i/j has a part SP 
of total probability transferred to 8 nearby lattice sites 
(e.g. at t = for tp the total probability is at a cer- 
tain lattice site, while for ip this site contains the prob- 
ability 1 — SP and nearby 8 sites have equal probabili- 
ties dP/8). For convenience we fit the decay of fidelity 
normalized by its maximal value fmax by the exponen- 



7 



tial decay f{t)/ fmax = exp(— Fi) with a certain decay 
rate F. The value of F is fixed by the condition that 
curves for various values of 5P are approximately super- 
imposed on one scaling curve as it is shown for example 
in Fig. [10] for W ~ % and /3 = 1. The same procedure 
was done for other values of disorder W . The result- 
ing dependence of F on 5P and the average IPR at 
large times at /3 = is shown in Fig. [TT] The data can 
be described by the dependence F (^-P/Co)^^^- We 
interpret this in the following way: the perturbationiJP 
spreads over states and gives a modification of the 
nonlinear potential ~ ((5P/^o)^^^ that determines 

a typical transition frequency to other states leading to 
F ~ /?(5|V'P l3{6P/ioY''^. A more detailed check of the 
functional dependence requires larger variation of that 
can be done in future studies. 
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FIG. 11: (Color online) Dependence of the fidelity decay rate 
r on rescaled perturbation 5P/^o for W — 8 (triangles), 10 
(squares), 15 (circles). The straight line shows the algebraic 
dependence given by fit: P = {5P/£,o)^ with — 0.486. 



the growth of the decay rate F becomes saturated and 
F reaches its saturated value F^ (see Fig. [121). The de- 
pendence of Tg on ^0 is shown in Fig. [T31 Except the 
strongly localized case at W = 15 the dependence is sat- 
isfactorily described by Tg ~ /3/Co- This corresponds 
to the situation when the perturbation of nonlinear field 
is rather strong and the decay of fidelity is given by a 
typical nonlinear frequency shift Sto ~ /Sji/'nP ~ 
Of course, this relation is valid on relatively short time 
scales used for investigation of fidelity decay (Figs. [TUl[T51) 
when IV'nP ^0- On a larger time scales the grows of 
^{t) with time should be taken into account. 




10 100 



^0 



FIG. 13: (Color online) Dependence of saturated decay rate 
Fs (like in Fig. ^ on the IPR for W = 15, 12, 10, 9, 8 
(from left to right). The slope of the straight line fit is 
-0.831 ±0.0095. 

The nonlinear fidelity decay gives new additional char- 
acteristics of nonlinear field evolution on moderate time 
scales. 



E 




0.5 1 1.5 2 

Ft 



FIG. 12: (Color online) Rescaled fidelity decay as a function 
of time, for VF = 8, /3 = 1 and values of 5P from 0.005 (top 
curve) to 0.7 (bottom curve); averaging is done over Nd = 12 
disorder realisations. At large 5P the behavior of /(t)//max 
becomes independent of SP. The straight line shows the decay 
with the saturated value of Fs'. f{t)/fmax ~ exp(— Fst). 

It is interesting to note that with the increase of 5P 



V. CONCLUSIONS 

The results of extensive numerical simulations pre- 
sented above show that at moderate nonlinearity in dis- 
ordered lattices with localized linear cigcnmodcs in ID 
and 2D there is an algebraic spreading over the lattice 
with the number of populated sites growing as An cx 
(see Eq. ([3])). This spreading continues up to enormously 
large times t = 10^ while for the linear problem the lo- 
calization takes place on a time scale tioc ~ 10 or 100. 
This result is in a satisfactory agreement with the pre- 
vious studies [2l|, [g^ . The numerical data are obtained 
on extremely large time scales (up to t = 10^ in dimen- 
sionless units) that are by 4 (Figs. to 8 (Figs. 
orders of magnitude larger than the time scale of Ander- 
son localization (down to tioc = 10). This indicates that 
the numerical results demonstrate the real asymptotic 
regime of algebraic growth. Even if the numerical simu- 
lations do not allow to make rigorous conclusions about 
the asymptotic spreading at infinite times we think that 
the enormous difference between tioc scale and the com- 
putational times reached in our numerical simulations 
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favors the conclusion that the presented numerical data 
give the real asymptotic behavior. The theoretical ex- 
ponent of spreading v = djiZd + 2) given by Eq. ([3]) 
is in a good agreement with the numerical data for the 
KNR model in ID ^ and the 2D Anderson model (P) 
(see Figs. I4l5l7l8p . For the ID Anderson model ^ the 
numerical value of the exponent v shows about 20-30% 
deviation from the theoretical value. We think that such 
a deviation should be attributed to a specific property 
of the ID Anderson model which has no diffusive regime 
showing a direct transition from a ballistic dynamics to 
localization (see e.g. [131 )• The data obtained for a weak 
nonlinearity with /? < 0.03 show no spreading up to times 
t = 10^10^ (see Fig. H Ref. [H, Fig. O Fig. [5] respec- 
tively). A small spreading seeing in the KNR model at 
enormously large time t = 10^ may indicate that a very 
slow (logarithmic ?) spreading in time is not completely 
excluded and processes like the Arnold diffusion [ij] may 
be present. However, this spreading, even if present, is 
so slow that in global the presented numerical data can 
be considered as a confirmation of the theoretical expec- 
tation according to which for a typical initial state the 



localization is preserved at /3 < /3c- Our data indicate 
that /3c ~ 1/30. Indeed, the spreading behavior is quali- 
tatively different for /? ^ 1 and /? ^ 1/30. 

It is possible that such type of slow probability and 
energy spreading over disordered lattices may play an 
important role in complex molecules giving more rapid 
propagation of probability and energy along molecular 
chains compared to a simple diffusion produced by noise. 
It would be interesting to observe the nonlinear destruc- 
tion of localization for BEG in disordered potential or 
for nonlinear waves in photonic lattices but this is rather 
hard task since very long observation times are required 
for that. 

We thank A.S.Pikovsky for stimulating discussions, 
one of us (DLS) thanks the participants of the NLSE 
workshop at the Lewiner Institute at the Technion for 
useful discussions. 

Note added: after the submission of this paper there 
appeared the preprint (44j where the same nonlinear ID 
Anderson model is investigated numerically and analyti- 
cally. 
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